clc
clear all

n = 500;
% m = 30; % order of cyclic group
p = 0.8;   % observation probability
q = 0.8; % probability of not being corrupted by (multiplicative) outlier noise, q=1 means no outlier noise

Number_of_m = 6;   % number of points
Array_m = 28 + 2* (1:Number_of_m)';   


Number_of_Sim = 5;   % number of repeats

% error 
MCR_Spectral = zeros(Number_of_Sim, Number_of_m);
MCR_PPM = zeros(Number_of_Sim, Number_of_m);
MCR_GPM = zeros(Number_of_Sim, Number_of_m);
MCR_Spectral0 = zeros(Number_of_Sim, Number_of_m);
MCR_GPM0 = zeros(Number_of_Sim, Number_of_m);

% running time
time_Spectral = zeros(Number_of_Sim, Number_of_m);
time_PPM = zeros(Number_of_Sim, Number_of_m);
time_GPM = zeros(Number_of_Sim, Number_of_m);
time_Spectral0 = zeros(Number_of_Sim, Number_of_m);
time_GPM0 = zeros(Number_of_Sim, Number_of_m);

for index_m = 1:Number_of_m
    m = Array_m(index_m);    
    
    for index_sim = 1:Number_of_Sim
        fprintf('index_n: %d of %d, index_sim: %d of %d\n', index_m, Number_of_m, index_sim, Number_of_Sim);
        
        %% Generating instance...
        % Set parameters
        param.n             = n;     % number of variables
        param.m             = m;      % number of states 
        param.p_obs         = p;       % observation ratio
        param.maxIter       = 50;      % maximum number of projected power iterations
        param.mu0           = 10;      % scale factor
        param.pi0           = q;     % non-corruption rate
        param.round_flag    = 'no';    % indicate whether mu0 is infinity
        param.max_poweriter = 20;      % maximum number of power iterations
        param.eps_power     = 0.02;    % tolerance level of power methods

        % Noise density for the random corruption model
        P0          = (1-param.pi0)* ones( param.m, 1) / param.m; 
        P0(1)       = P0(1) + param.pi0;
        param.P0    = P0;        
        % display(sprintf('n=%d, m=%d, corruption_rate=%f', param.n, param.m, 1-param.pi0));

        % Generate a synthetic example
        [Problem, C, S_true]     = ProblemGenerator(param);
        

        %% Projected power method
        tic;
        L           = Problem.L;        % Input log-likelihood matrix
        Lambda      = eigs( L, 2);      % Compute 2nd largest eigenvalue of L
        L_opt       = @(I) L  * I;      % Matrix-vector multiplication operator
        proj        = @(z) ProjectOperator(z, param.n, param.m, param.mu0/Lambda(2)); % Projection operator
        [z, relErrors, z0, t_Spectral0]  = noncvxAlign(param, Problem.X_gt, L_opt, proj);  % z is the output and relErrors is the L2 error
        S_Spectral0 = z2S(z0,m);
        t_PPM = toc;
        S_PPM = z2S(z, m);
        
        %% Spectral Estimator
        [t_Spectral, S_Spectral] = Spectral_Cyclic(C, m);
        
        %% GPM
        [t_GPM, S_GPM] = GPM_Cyclic( C, S_Spectral, m );
        
        %% GPM0
        [t_GPM0, S_GPM0] = GPM_Cyclic( C, S_Spectral0, m );
        
        %% Update error
        MCR_Spectral(index_sim, index_m) = MCR(S_true, S_Spectral, m);
        MCR_Spectral0(index_sim, index_m) = MCR(S_true, S_Spectral0, m);
        MCR_GPM(index_sim, index_m) = MCR(S_true, S_GPM, m);
        MCR_GPM0(index_sim, index_m) = MCR(S_true, S_GPM0, m);
        MCR_PPM(index_sim, index_m) = MCR(S_true, S_PPM, m);
        
        %% Update time
        time_Spectral(index_sim, index_m) = t_Spectral;
        time_Spectral0(index_sim, index_m) = t_Spectral0;
        time_GPM(index_sim, index_m) = t_GPM + t_Spectral;
        time_GPM0(index_sim, index_m) = t_GPM0 + t_Spectral0;
        time_PPM(index_sim, index_m) =  t_PPM;
    end
end

% 
fig = figure;
set(gcf,'units','points','position',[10,10,800,400])
% title('d=5, observe p_1 = 0.1, noise p_2=0.7')

subplot(1,2,1);
set(gca,'fontsize',15)
plot(Array_m, 1 - sum(MCR_Spectral)/Number_of_Sim,'-kp','LineWidth',2,'MarkerSize',10);
hold on
plot(Array_m, 1 - sum(MCR_GPM)/Number_of_Sim,'-ro','LineWidth',2,'MarkerSize',10);
plot(Array_m, 1 - sum(MCR_PPM)/Number_of_Sim,'-g*','LineWidth',2,'MarkerSize',10);
plot(Array_m, 1 - sum(MCR_Spectral0)/Number_of_Sim,'-mx','LineWidth',2,'MarkerSize',10);
plot(Array_m, 1 - sum(MCR_GPM0)/Number_of_Sim,'-b+','LineWidth',2,'MarkerSize',10);
hold off
legend('Spectral','GPM','PPM','Spectral_{cc}','GPM_{cc}','Location','southwest')
xlabel('$m$','Interpreter','latex')
ylabel('\textbf{Recovery Rate}','Interpreter','latex')


subplot(1,2,2);
set(gca,'fontsize',15)
plot(Array_m,sum(time_Spectral)/(Number_of_Sim),'-kp','LineWidth',2,'MarkerSize',10);
hold on
plot(Array_m,sum(time_GPM)/(Number_of_Sim),'-ro','LineWidth',2,'MarkerSize',10);
plot(Array_m,sum(time_PPM)/(Number_of_Sim),'-g*','LineWidth',2,'MarkerSize',10);
plot(Array_m,sum(time_Spectral0)/(Number_of_Sim),'-mx','LineWidth',2,'MarkerSize',10);
plot(Array_m,sum(time_GPM0)/(Number_of_Sim),'-b+','LineWidth',2,'MarkerSize',10);
hold off
legend('Spectral','GPM','PPM','Spectral_{cc}','GPM_{cc}','Location','northwest')
xlabel('$m$','Interpreter','latex')
ylabel('\textbf{Time (s)}','Interpreter','latex')

% sgtitle(['$n$ = ', num2str(n), ', $p$ = ', num2str(p), ', $q$ = ', num2str(q)],'Interpreter','latex');

% % Save
% folder_name = ['n=', num2str(n),'_m=',num2str(Array_m(1)),'-',num2str(Array_m(Number_of_m)), '_p=', strrep(num2str(p),'.',''), '_q=', strrep(num2str(q),'.','')];
% mkdir('C:\Users\AMA\Dropbox\Academics and Research\Research Projects\Working\Synchronization over Orthogonal Group\experiment2020\Cyclic\Results',folder_name);
% save(['Results\',num2str(folder_name),'\MCR_Spec.mat'], 'MCR_Spectral');
% save(['Results\',num2str(folder_name),'\MCR_GPM.mat'], 'MCR_GPM');
% save(['Results\',num2str(folder_name),'\MCR_PPM.mat'], 'MCR_PPM');
% save(['Results\',num2str(folder_name),'\time_Spec.mat'], 'time_Spectral');
% save(['Results\',num2str(folder_name),'\time_GPM.mat'], 'time_GPM');
% save(['Results\',num2str(folder_name),'\time_PPM.mat'], 'time_PPM');
% savefig(['Results\',num2str(folder_name),'\',num2str(folder_name),'.fig']);